Moment ratios for the pair contact process with diffusion 
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Abstract 

We study the continuous absorbing-state phase transition in the one-dimensional pair contact 
process with diffusion (PCPD). In previous studies [Dickman and de Menezes, Phys. Rev. E, 
66 045101 (R) (2002)], the critical point moment ratios of the order parameter showed anomalous 
behavior, growing with system size rather than taking universal values, as expected. Using the 
quasistationary simulation method we determine the moments of the order parameter up to fourth 
order at the critical point, in systems of up to 40960 sites. Due to strong finite-size effects, 
the ratios converge only for large system sizes. Moment ratios and associated order-parameter 
histograms are compared with those of directed percolation. We also report an improved estimate 
[p c = 0.077092(1)] for the nondiffusive pair contact process. 
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I. INTRODUCTION 

The pair contact process (PCP) [ill is a nonequilibrium stochastic model that exhibits a 
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phase transition to an absorbing state |2J, y, LI ■ Several studies have established that the 
PCP belongs to the robust directed percolation (DP) universality class j^. In contrast with 
the contact process, which has the vacuum as its unique absorbing state, the PCP exhibits 
an infinite number of absorbing configurations, since both creation and annihilation require 
a nearest-neighbor pair of particles. Allowing particles in the PCP to hop on the lattice, 
one obtains the so-called pair contact process with diffusion (PCPD). Here there are only 
two absorbing states: the vacuum, and the subspace with only a single particle. 

While a version of the pair contact process with diffusion was proposed by Grassberger in 
1982 6], current interest in the problem follows its rediscovery by Howard and Taiiber p| who 
suggested that its Langevin description involves complex noise. The first numerical studies 
of the PCPD suggested that critical the PCPD would fall in the "parity conserving" 
(PC) class, but increasing computational effort revealed that this was not the case, and 
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that the critical behavior of the model was masked by huge correct ions |9J, llfj. Since then, 
diverse scenarios, some of them contradictory, have been proposed in order to clarify this 
question. At present there are two principal schools of thought: in one, the PCPD belongs 
to a novel universality class distinct from DP, with a unique set of critical exponents, or 
possibly continuously varying exponents due to a marginal perturbation Q,Q,Q. The 
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ds that the PCPD should be attracted to a DP fixed-point after a huge 



15 1 . A recent review on these and other scenarios can be found in |l6| . 



In this work we study the PCPD via quasistationary (QS) simulations 13], focussing on the 
order parameter and its moments. 

The balance of this paper is organized as follows: In the next section we review the 
definition of the model and detail our simulation method. In Sec. Ill we present our results, 
and Section IV is devoted to discussion and conclusions. 



II. MODEL AND SIMULATION METHOD 



The PCP is defined on a lattice, with each site either occupied or vacant. All changes of 
configuration involve a pair of particles occupying nearest-neighbor sites, called a pair in what 



follows. A pair annihilates itself at a rate of p, and with rate 1—p creates a new particle at a 
randomly chosen site neighboring the pair, if this site is vacant. In the PCPD, in addition to 
the creation and annihilation processes already mentioned, each particle attempts to hop, at 
rate D, to a randomly chosen nearest- neighbor (NN) site; the move is accepted if the target 
site is vacant. The PCPD exhibits a continuous phase transition to the absorbing state, 
at a critical annihilation rate p c (D). Several variants of the model, differing in how each 
process (creation, annihilation or diffusion) is selected, have been studied in the literature 
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13. Here we use the implementation of Ref. The model is defined on a ring of L 

sites. We maintain a list of pairs to improve efficiency. At each step of the evolution we first 
choose between diffusion (with probability D) and reaction (with probability 1 — D). In a 
reaction step, we select a pair at random, and choose between creation and annihilation with 
probabilities p and 1—p, respectively. If we choose diffusion, a particle is selected at random; 
it attempts to hop to one of its neighbor sites (if the latter is vacant). The time increment 
associated with each step is At = l/(N pair + DN part ), where N pair and N part are the number 
of pairs and particles in the lattice at time t. Thus each lattice site is effectively visited 
once (on average) per time unit. The most obvious definition of the order parameter in the 
PCPD is the pair density p2, the number of pairs per site. Since creation (and destruction) 
of particles requires pairs, one might expect the particle density p\ to scale in a similar 
manner. 

In the studies reported here we sample the quasistationary (QS) distribution of the pro- 
cess, (that is, conditioned on survival), which is very useful in the study of processes with 
an absorbing state. (In fact, conventional simulations of "stationary" properties of lattice 
models with an absorbing state actually study the quasistationary regime, given that the 
only true stationary state for a finite system is the absorbing one). We employ a recently 
devised simulation method that yields quasistationary properties directly [l7|. This is done 
by maintaining, and gradually updating, a set of configurations visited during the evolu- 
tion; when a transition to the absorbing state is imminent the system is instead placed in 
one of the saved configurations. Otherwise the evolution is exactly that of a conventional 
simulation. 

The above scheme was shown ^1 to yield precise results, in accord with the exact QS 
distribution for the contact process on a complete graph, and with conventional simulations 
of the same model on a ring jl7j |. The scheme has also been shown to yield results that 



agree, to within uncertainty, with the corresponding results of conventional simulations for 
a sandpile model . The advantage of the method is that a realization of the process can 
be run to arbitrarily long times. Thus, whereas in conventional simulations a large number 
of realizations must be performed to have a decent sampling of the (quasi)stationary state, 
here every realization provides useful information, once the initial transient has relaxed. 
This leads to an order of magnitude improvement in efficiency, in the critical region. For 



further details on the method see 
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The PCPD dynamics is characterized by (1) single-particle diffusion and (2) reactions 
involving a pair, motivating us to ignore the "purely diffusive" subspace. This means that 
we modify slightly the dynamics of the model, restricting it to the subspace with at least 
one pair, which we call the reactive subspace. The motivations for studying the modified 
process are: (1) In |llj, excluding the subspace without pairs yielded better behaved moment 
ratios; (2) eliminating the large fraction of time expended on the nonreactive intervals yields 
a major improvement in efficiency We may further justify our choice by noting that any 
scaling properties necessarily involve reactions, that is, the existence of pairs. 

We impose this restriction as follows. The initial configuration (all sites occupied) has a 
large number of pairs. Applying the QS simulation method, we accumulate a list of config- 
urations during the evolution. Then, whenever a visit to a configuration (absorbing or not) 
without pairs is imminent, the system is instead placed in a configuration selected randomly 
from the list. In other words, the QS method, previously used to sample quasistationary 
properties, is used here to restrict the dynamics to the subspace of interest. A subtle differ- 
ence between the two applications of the method is that, while in usual QS simulations 17 1 
the method provides a just sampling of properties conditioned on survival, in the present 
case we eliminate certain nonabsorbing configurations (and the histories involving them) as 
well. For the reasons given above, we do not expect this to color our results for large system 
sizes. 



III. SIMULATION RESULTS 

We performed extensive simulations of the PCPD on rings of L = 1280, 2560, ...40960 
sites, using the QS method restricted to the reactive subspace. Each realization of the 
process is initialized with all sites occupied, and runs for 10 9 time steps. Averages are taken 



in the QS regime, after discarding an initial transient which depends on the system size. 
In practice we accumulate histograms of the time during which the system has exactly 1, 
2,...n,... pairs, and similarly for particles. The histograms are used to evaluate moments; we 
denote by m J;1 the j-th moment of the particle number probability distribution; 771,- 2 denotes 
the corresponding moment of the pair number distribution. (Thus the order parameter p 2 
could also be denoted mi ;2 .) The QS lifetime r is taken as the mean time between attempts 
to leave the reactive subspace. 

The number of saved configurations ranges from 10000, for L = 1280, to 500 for L = 
40960. Values of p rep range from 10~ 4 to 2 x 10~ 6 . The results of QS simulations were found 
to agree with results of conventional simulations h|, for system sizes L = 80, 160. ..1280. 
We study three diffusion rates, D = 0.1, 0.5 and 0.85, as in [111 ]. For comparison, we also 
study the (nondiffusive) PCP, whose scaling properties are known to be those of directed 
percolation. The QS results for the D = case confirm DP behavior in the nondiffusive 
PCP, as can be verified from the values of the exponents and moment ratios listed on Tables 
I and II. Further, we improved the estimate for the critical point of the PCP, obtaining 
p c = 0.077092(1). (This is consistent with the previous best estimate, p c = 0.077090(5), of 
Ref. H-) 

The first step in analyzing our results is to determine, for each D value studied, the 
critical annihilation probability p(D). We use the following criteria for criticality: power- 
law dependence of (1) p\ and p 2 and (2) r on system size L (i.e., the usual finite-size scaling 
relations p ~ L - / 3 /^ and r ~ L z ); (3) constancy of the moment ratio r 2 n ; 2 = Tn^.^/m^.^ 
with system size. The three criteria were found to be mutually consistent within the error 
margins. Figure 1 shows the data for p 2 for D = 0.5; the QS order parameter for various 
p values near p c is plotted versus L on log scales. The data for the four largest sizes are 
well fit by a straight line of slope —/3/u± = —0.385. Plotting L l3 ^ v± p 2 (see inset), allows 
us to eliminate as off-critical p values for which the plot shows a significant curvature, 
leading to the estimate p c = 0.120353(2). A similar analysis of the particle density p\ yields 
p c = 0.120357(5) while that for the lifetime r gives 0.120352(3), leading to the overall best 
estimate p c 0. 120354(3) for D = 0.5. 

Analysis the moment ratio is also useful in setting limits on p c : As can be seen in Fig. 2, 
this quantity appears to grow with system size for p < p c , and vice-versa. For D = 0.85 
for example, we find p c = 0.12992(1) from analysis of p 2 > Pc — 0.12993(1) from analysis of 



r, and p c = 0.129925(8) from analysis of the moment ratio r 2 u-2- (For D = 0.5, on the 
other hand, the moment ratio data do not yield an estimate of precision comparable to that 
furnished by pi, p 2 and r.) 

As a check on our procedure for determining p c , we also performed, for D = 0.85 initial 
decay stu.es HQ. lo these studies .eo^p—; p 2 is followed function of 
time, starting, as in the other studies, with all sites occupied. (In this case the system does 
not leave the reactive subspace on the time scale of the simulation so the QS procedure is not 
needed.) Here the expected behavior is p 2 ~ t~ e . Using deviations from the power law to 
identify off-critical values, a study of systems of 10 6 sites, to a maximum time of 10 8 , yields 
p c = 0.129915(15), fully consistent with the QS results. Analysis of the data for t > 5 x 10 5 
furnishes 9 = 0.19(1), consistent with previously reported results jl^l l^. 

Our present estimates for p c are slightly lower than those reported in [ijj, a difference 
of less than 0.1%. These differences highlight the strong finite-size corrections affecting the 
PCPD (in [ljj system sizes range from 80 to 1280), and appear in other absorbing-state 
problems, such as the restricted sandpile model [lsj ]. 

Using our results for the three largest system sizes (L = 10240, 20480 and 40960) we 

obtain estimates for the critical exponent ratios (3/v±_ and z = v\\/v±- These are reported in 

Table I. The chief contribution to the uncertainties in the exponent ratios comes from the 

uncertainty in p c . We note that although the exponent z appears to depend systematically 

on diffusion rate D, our results are in fact consistent with z = 2 in all cases. By contrast, 

the value of fi/u± for D = 0.1 is significantly different than that found for D = 0.5 and 0.85. 

The latter value {(3/v±_ = 0.385(11)) does not yield a good fit to the data for D = 0.1 even 

if a logarithmic correction term is introduced in the fitting function. 

In nonequilibrium statistical physics, obtaining values for moment ratios has proved an 
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efficient method for identifying the universality class [19(. Here we analyze, in addition to 
r 2ii ; 2 mentioned above, the ratios 

r 3 i2-,2 = (1) 

m 2; 2mi ;2 

and 

^42;2 = — q— (2) 
m 2;2 

and the corresponding ratios for particles. Scaling arguments imply that such ratios such 
assume universal values at the critical point, as has been amply verified in equilibrium, and 



for models with an absorbing state such as the contact process and the PCP 



19j | . Analysis of 



such ratios is of particular interest in the present context, given the perplexing result of Ref. 



11] , namely, apparently unlimited growth in r 2 n ; 2 with increasing system size, suggesting 
that the PCPD does not follow the usual scaling behavior observed at absorbing-state phase 
transitions. The present large-scale simulations show that the PCPD moment ratios do in 
fact exhibit the expected behavior. 

Figure 3 shows the moment ratios for pairs and particles in the critical PCPD, with 
D = 0.5 and the PCP. On the basis of these data we estimate lini^oo T2\i-2 = 1.166(8) for 
the PCPD with D = 0.5. While the results for r 2 u-2 appear to converge for the system sizes 
studies, the data for the particle moment ratio r 2 u-i continue to grow with system size and 
cannot be extrapolated to the infinite-size limit. (We note in passing that the values for 
particles and pairs suggest that the two ratios become equal at a system size on the order of 
5 x 10 5 , for D = 0.5. Since universality implies similar scaling of particle and pair densities 
at the critical point, we may take this as an order of magnitude estimate for the system size 
at which corrections to scaling are finally superseded. The moment ratio data for D = 0.85 
are consistent with this estimate, while those for D = 0.1 suggest convergence only at even 
larger sizes, on the order of L ~ 10 7 .) Table II lists estimates of the the limiting (L — > oo) 
moment ratios for the PCPD and some other known universality classes. The ratios for 
D > 0.5 are rather similar to those for directed percolation. The cumulant ratio K±jK\ 
(here Kj denotes the j-th cumulant of the order parameter probability distribution) is not 
in good accord with the DP value, although the estimates for he PCPD are rather imprecise. 

Another quantity that is useful for determining the universality class is the scaled order 
parameter histogram, defined as p*(n*) = np(n), where n denotes the number of pairs and 
n its mean value. (Here n* = n/n.) Figure 4 compares the scaled histograms for the PCPD 
with D = 0.1, 0.5 and 0.85. The latter two are quite similar, while the curve for low diffusion 
rate has a somewhat narrower peak. In the inset of Fig. 4 we compare the scaled histogram 
for the PCPD with D = 0.5 to that of the critical contact process and PCP. While the 
CP and PCP have virtually identical histograms, it is clear that even for the largest sizes 
studied the PCPD histogram is very different. (A striking difference between the PCPD 
on one hand and CP /PCP on the other is the behavior near n* = 0: in the latter case the 
probability increases linearly with n, while former exhibits a distinctly parabolic shape. The 
parabolic form is not an artefact of the QS simulation method, as it is already observed in 



Ref. using conventional simulation.) 

In the case of D = 0.85, we accumulated data for several values of the annihilation rate p 
in the vicinity of p c . This permits us to estimate the correlation length exponent u±, using 
the finite-size scaling relation 

m(A,L) oc J r m (L 1 ^ u± A), (3) 
where A = p — p c and T m is a scaling function. This leads to 

ocL 1 ^. (4) 



dm 



The data of Fig. 5 yield the estimate z/j_ = 1.09(1), confirming the results of [11 
IV. DISCUSSION 

We perform large-scale simulations of the pair contact process with diffusion, and find 
that ratios of various moments of the order parameter approach finite limiting values for 
large system sizes, allowing us to present the first numerical estimates of such ratios for 



this model. This resolves the apparent anomaly reported in [11 1. of a moment ratio growing 
without limit. 

Two basic questions prominent in recent discussions of the PCPD are: (1) Does the model 
exhibit continuously variable critical exponents, or a unique set, independent of diffusion 
rate D? (2) In the latter the exponents those of directed percolation? Despite the 

large computational effort, we are unable to answer these questions definitively. Our results 
nevertheless provide some clues. With regard to the first question, the results reported in 
Tables I and II show very good agreement for the two larger diffusion rates studied (D = 0.5 
and 0.85), but those for D = 0.1 are quit e different. Our values for the critical exponents 
are consistent with earlier studies ([ill ll^ for low D and [ill. for high D). 

For low diffusion rate (D — 0.1) we find evidence of stronger corrections to scaling, and/or 
a slower rate of convergence with increasing L, as noted above in the analysis of the moment 
ratios. Our results of z = 2.08(15) and /3/u± = 0.505(10) are consistent with previous results 



for low-diffusion rate 
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, 3- While tins argues in favor of nonuniversal behavior, we cannot 
rule out the possibility of a unique set of exponents, as suggested for example by Odor [22] , 
who asserted that, including a logarithmic correction, the same value for /3/v\\ holds for both 



high and low diffusion rates. (We note, however, that we were unable to fit all the data 
using the same value for [3/v±, even including a logarithmic correction.) In our opinion, it 
is still unclear if the observed difference in the exponents and moment ratios for D = 0.1, as 
compared with larger diffusion rates, is due to corrections to scaling. Such corrections would 
have to be exceptionally large to account for the observed differences, and have to exhibit a 
rather irregular behavior, given the close agreement in the results for D = 0.5 and 0.85. An 
equally if not more natural interpretation is that the critical exponents and moment ratios 
really do depend on D. 

If we accept, provisionally, that the PCPD critical properties are independent of D (with 
huge corrections to scaling for small diffusion rates), it is still not obvious that these prop- 
erties are characteristic of the DP universality class. The values for (3/v±_ and z are quite 
far from those of DP, even for larger diffusion rates, where corrections to scaling appear to 
be less severe. It is true that our estimates for u± (for D = 0.85), and for the moment ratios 
(for both D = 0.5 and 0.85) are close to the DP values, but the qualitatively different form 
of the order parameter histograms again prevents identification of the PCPD as belonging to 
the DP universality class, even for higher diffusion rates. Thus our results tend to support 
the conclusion of Ref. |2l|, that PCPD scaling is distinct from that of DP. 

Large-scale studies of time-dependent behavior (2QJ suggested a value for the initial decay 
critical exponent 9 different from that of DP. A subsequent study [3| led to the suggestion 
that apparent values of critical exponents vary as one increases the simulation time and 
system size, reflecting strong corrections to scaling, and that observed exponent values be 
interpreted as upper limits on the true values. The results of the quasistationary simulations 
are not affected by finite-time corrections, although finite size effects are still present. 

In summary, we have applied the quasistationary simulation method to the pair contact 
process with diffusion, restricted to the subspace with at least one pair. Our results indicate 
that the anomalous behavior of the critical order-parameter moment ratios does not persist 
for large systems. Restricting the dynamics to the reactive sector, these ratios appear to 
converge to finite values, as expected. Taken at face value our results imply a variation 
of scaling properties with diffusion rate, but the opposite interpretation is tenable if one 
invokes strong corrections to scaling for small D. Regardless of whether or not the critical 
exponents vary with D, it seems premature to conclude that the PCPD will eventually cross 
over to the DP class. 
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Table I. Critical exponent values for the PCP, PCPD and DP. DP values from Ref. [2 



D V c 



f3/v± z 

1.584(7) 
2.08(15) 
2.04(5) 
1.88(12) 
1.5807(1) 



(PCP) 0.077092(1) 0.2519(3) 
0.1 0.106405(15) 0.505(10) 
0.5 0.120354(3) 0.385(11) 
0.85 0.129925(8) 0.386(5) 
DP - 0.25208(5) 



Table II. Critical moment ratio values for the PCP, PCPD and the DP, parity-conserving (PC), 

and conserved-DP (C-DP) universality classes. 



D 


?"211;2 ?"312;2 


f 422;2 


K±jK\ Reference 


(PCP) 1.1738(2) 1.303(3) 


1.558(2) 


-0.493(3) [12] 


0.1 


1.140(15) 1.27(2) 


1.55(3) 


0.1(2) this work 


0.5 


1.166(8) 1.310(15) 1.61(2) 


0.0(1) this work 


0.85 


1.170(6) 1.31(1) 


1.61(4) 


-0.1(1) this work 


DP 


1.1736(2) 1.301(3) 


1.554(2) 


-0.505(3) [19] 


PC 


1.3340(4) 




[11] 


C-DP 


1.142(8) 







FIGURE CAPTIONS 



FIG. 1. Quasistationary order parameter versus system size for p = 0.120345, p = 0.120350, 
p = 0.120355 and p = 0.120360, from top to bottom. D = 0.85. Inset: lnL^p versus InL 
for the same values of p. 

FIG. 2. Quasistationary moment ratio r 2 n ;2 versus system size for p = 0.12990, p = 0.12993, 
and p = 0.12995, from top to bottom D = 0.85. 

FIG. 3. Moment ratios r 211;2 (pairs) and r 2 n ; i (particles) versus InL for p = 0.120354 and 
D = 0.5. Dashed line: moment ratio r 211;2 for the critical PCP. 

FIG. 4. Scaled histograms for the critical PCPD with D — 0.1, 0.5, and 0.85 (from top to 
bottom). Inset: scaled histograms for the critical PCPD (D = 0.5, L = 10240 and 20480), 
the critical CP (same sizes) and the critical PCP (L = 10240). 

FIG. 5. Quasistationary moment ratio r 2 n ; 2 versus p for system sizes L = 1280, L = 
2560. ..L = 20480 and D = 0.85. Inset \ndm/dp versus InL at criticality. 
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